Set wd and load packages.
Intensity data is handled here.
The desired spatial extent can be set via subsetting: [extent_object].
dates <- c(ymd("2018-01-03"), ymd("2018-01-15"), ymd("2018-01-27"), ymd("2018-02-08"), ymd("2018-02-20"), ymd("2018-03-04"), ymd("2018-03-16"), ymd("2018-03-28"), ymd("2018-04-09"), ymd("2018-04-21"), ymd("2018-05-03"), ymd("2018-05-15"), ymd("2018-05-27"), ymd("2018-06-08"), ymd("2018-06-20"), ymd("2018-07-02"), ymd("2018-07-14"), ymd("2018-07-26"), ymd("2018-08-07"), ymd("2018-08-19"), ymd("2018-08-31"), ymd("2018-09-12"), ymd("2018-09-24"), ymd("2018-10-06"), ymd("2018-10-18"), ymd("2018-10-30"), ymd("2018-11-11"), ymd("2018-11-23"), ymd("2018-12-05"), ymd("2018-12-17"), ymd("2018-12-29"))
int_vv <- st_set_dimensions(int_vv, 3, values = dates, names = "time")
int_vh <- st_set_dimensions(int_vh, 3, values = dates, names = "time")
Needs to be freshly warped for some reason, can’t be saved (write_stars) and then loaded. But: rain_warped.tif is a file with already aggregated rain data (31 time steps instead of 365).
Then combine all three. Give along = "bands" so that all attributes are treated as bands.
This is more of a design choice.
comb_split <- split(comb, "bands")
names(comb_split) <- c("VV", "VH", "prec")
comb_split
## stars object with 3 dimensions and 3 attributes
## attribute(s):
## VV VH prec
## Min. :-20.54 Min. :-24.74 Min. : 0.00
## 1st Qu.:-15.68 1st Qu.:-20.72 1st Qu.: 0.00
## Median :-14.50 Median :-19.72 Median : 0.20
## Mean :-14.17 Mean :-19.59 Mean : 2.59
## 3rd Qu.:-12.66 3rd Qu.:-18.45 3rd Qu.: 3.40
## Max. : -5.13 Max. :-10.62 Max. :17.60
## NA's :162998 NA's :162998 NA's :162998
## dimension(s):
## from to offset delta refsys point values x/y
## x 4745 4865 715082 10 ETRS89 / UTM zone 32N FALSE NULL [x]
## y 1083 1202 5864854 -10 ETRS89 / UTM zone 32N FALSE NULL [y]
## time 1 31 2018-01-03 12 days Date NA NULL
# get one pixel ts (at 10/10)
pixel.df <- as.data.frame(comb_split[1,100,100,])
pixel.df <- cbind(pixel.df, comb_split[2,100,100])
pixel.df <- cbind(pixel.df, comb_split[3,100,100])
pixel.df <- pixel.df[,c(3,4,8,12)]
# calculate additional measures
names(pixel.df) <- c("time", "VV", "VH", "prec")
pixel.df$sum <- pixel.df$VV + pixel.df$VH
pixel.df$diff <- pixel.df$VV - pixel.df$VH
pixel.df$ratio <- (pixel.df$VV / pixel.df$VH)
# plot all
ggplot() +
geom_bar(aes(x=pixel.df[,1], y=pixel.df$prec), stat='identity') +
geom_line(aes(x=pixel.df[,1], y=pixel.df$VV, color="VV")) +
geom_line(aes(x=pixel.df[,1], y=pixel.df$VH, color="VH")) +
geom_line(aes(x=pixel.df[,1], y=pixel.df$sum, color="sum")) +
geom_line(aes(x=pixel.df[,1], y=pixel.df$diff, color="diff")) +
geom_line(aes(x=pixel.df[,1], y=pixel.df$ratio, color="ratio")) +
scale_color_manual(name="Polarization", values=c(VV="red", VH="blue", sum="lightgreen", diff="forestgreen", ratio="yellow"))
# take a look at the ratio
ggplot() +
geom_line(aes(x=pixel.df[,1], y=pixel.df$ratio))
A decrease in VV together with not so much decrease in VH might point to flooding with intact vegetation. More information on the plant appearance is needed.
Giving the opportunity of constructing very different results. Depending on whether neighboring pixels or random pixels are selected, similarities and differences are seen.
# make first attribute, pixel 10/10 into dataframe template
pix.df <- as.data.frame(comb_split[1,100,100,])[,3:4]
# make data frame for all pixels
for (i in 1:120) {
for (j in 1:119) {
if (is.na(comb_split[1,i,j,1][[1]])) {
}
else {
pix.df <- cbind(pix.df, as.vector(comb_split[1,i,j,][[1]]))
}
}
}
# delete template data frame column
pix.df <- pix.df[,c(1,3:9264)]
# assign numbers as names
names(pix.df) <- c("time", seq(1,9262,1))
Here we have the possibility to select which pixels are plotted. For example each 37th pixel in the dataframe:
# initiate plot
p <- ggplot(data=pix.df, aes(x=pix.df$time)) + ylab("VV")
# make sequence of pixel positions to plot
ind <- seq(2,9263, 67)
# loop
for (i in ind) { p <- p + geom_line(aes_string(y = pix.df[,i]), alpha = 0.2) }
p
or the first X columns, meaning a profile through the image:
p <- ggplot(data=pix.df, aes(x=pix.df$time)) + ylab("VV")
ind <- seq(2,60, 1)
for (i in ind) { p <- p + geom_line(aes_string(y = pix.df[,i]), alpha=0.3) }
p
Instead of making a dataframe containing VV backscatter for each timestep, a VV/VH ratio is calculated per pixel, per timestep.
VV and VH seem to grow more similar during the course of the year. Surprising unison in first less and then more similarity is exhibited around the start of april.
NDVI_0204 <- read_stars("./S2_data/NDVI_02042018_clip.tif")
MDWI_3_8 <- read_stars("./S2_data/MDWI_3_8_02042018_clip.tif")
MDWI_8_12 <- read_stars("./S2_data/MDWI_8_12_02042018_clip.tif")
shape_NDVI <- st_transform(shape, crs=st_crs(NDVI_0204))
NDVI_0204_clip <- st_as_stars(NDVI_0204[shape_NDVI[124,]])
MDWI_3_8 <- st_as_stars(MDWI_3_8[shape_NDVI[124,]])
MDWI_8_12 <- st_as_stars(MDWI_8_12[shape_NDVI[124,]])
# plot(comb_split[1,,,9], main="VV 09.04.")
# plot(NDVI_0204_clip, main="NDVI 02.04.")
# plot(comb_split[2,,,9], main="VH 09.04.")
plot(comb_split[1,,,9], main="VV 09.04.")
plot(comb_split[2,,,9], main="VH 09.04.")
# NDVI from 01.02.2018, S1 from 27.01.
NDVI_0102 <- read_stars("./S2_data/NDVI_01022018_clip.tif")
NDVI_0102_clip <- st_as_stars(NDVI_0102[shape_NDVI[124,]])
plot(NDVI_0102_clip, main="NDVI 01.02.")
plot(comb_split[1,,,3], main="VV 27.01.")
plot(comb_split[2,,,3], main="VH 27.01.")
dates[16:18]
## [1] "2018-07-02" "2018-07-14" "2018-07-26"
Sentinel Playground screenshot (NDWI - plant moisture) in clarification of the dark spot in the middle left part of the polygon. NDVI shows slightly less Vegetation. Visual shows a slightly brighter patch. -> not water.